Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 2 June 2009 (MN BTeX style file v2.2) 



Discovery of Drifting High-frequency QPOs in Global 
Simulations of Magnetic Boundary Layers 

M. M. Romano va*, A. K. Kulkarnif 

Dept. of Astronomy, Cornell University, Ithaca, NY 14853 
2 June 2009 

ABSTRACT 

We report on the numerical discovery of quasi-periodic oscillations (QPOs) associated 
with accretion through a non-axisymmetric magnetic boundary layer in the unstable 
regime, when two ordered equatorial streams form and rotate synchronously at ap- 
proximately the angular velocity of the inner disk The streams hit the star's surface 
producing hot spots. Rotation of the spots leads to high-frequency QPOs. We per- 
formed a number of simulation runs for different magnetospheric sizes from small to 
tiny, and observed a definite correlation between the inner disk radius and the QPO 
frequency: the frequency is higher when the magnetosphere is smaller In the stable 
regime a small magnetosphere forms and accretion through the usual funnel streams 
is observed, and the frequency of the star is expected to dominate the lightcurve. 
We performed exploratory investigations of the case in which the magnetosphere be- 
comes negligibly small and the disk interacts with the star through an equatorial belt. 
We also performed investigation of somewhat larger magnetospheres where one or 
two ordered tongues may dominate over other chaotic tongues. In application to mil- 
lisecond pulsars we obtain QPO frequencies in the range of 350 Hz to 990 Hz for one 
spot. The frequency associated with rotation of one spot may dominate if spots are not 
identical or antipodal. If the spots are similar and antipodal then the frequencies are 
twice as high. We show that variation of the accretion rate leads to drift of the QPO 
peak. 

Key words: accretion, dipole — plasmas — magnetic fields — stars 



1 INTRODUCTION 

In disk-accreting magnetized stars, disk-magnetosphere in- 
teraction and channeling of matter to the magnetic poles 
may determine the majority of observational features of such 
stars (e.g., |Ghosh & Lamb|1979l|K6nigl|1991|l. The se include 
young solar-tj^e stars (e.g., |Bouvier et a l. ||2007|l an d very 
old stars such as magnetized white dwarfs ( Warner 1995 ) or 
accreting millisecond pulsars ( van der Klis 2000) . In many 
cases no clear period is observed, and high-frequency quasi- 
periodic oscillations (QPOs) with a frequency corresponding 
to a fraction of the Keplerian frequency at the stellar surface 
are seen instead. It is often the case that a second, lower fre- 
quency peak appears and drifts in concert with the main fre- 
quency, forming two-peak QPOs. Such high-frequency oscilla- 
tions are observed in both accreting millisecond pulsars (e.g., 
[Wi jnands & van der Klis|19 98 ) and in dwarf novae, which are 
a sub-class of cataclysmic variables ( jWoudt & Warner|[2002[ 
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[Warner & Woudt l2006') . It has been suggested that the high- 
frequency QPOs may be associated with accretion to a star 
with a very weak magnetic field, when some sort of weakly 
magnetized equatorial belt forms (Paczyhski'1978'; Warner &J 
Woudt 2002). Miller et al. (1998) and Lamb & Miller (2001) 
suggested that the high-frequency QPOs in accreting millisec- 
ond pulsars may result from rotation of clumps in the inner 
disk around a weakly magnetized star. They suggested that 
some of the disk matter may penetrate through the magne- 
tosphere until stopped, e.g., by the radiation pressure. In this 
paper we investigate accretion to stars with small magneto- 
spheres, where the energy density of the disk matter domi- 
nates up to very small distances from the star, or even up to 
the stellar surface, while magnetospheres are small compared 
to the radius of the star. We call the region around the star 
with a small magnetosphere the Magnetic Boundary Layer 
(MBL), and consider different cases. We notice that as for 
large magnetospheres, matter may accrete either in the sta- 
ble (e.g.,|Romanova et al. 2003 , 2004, hereafter RUKL03 and 
RUKL04) or unstab le regime (Kulkarni & Romanova 2008a b| 
hereafter KR08a,b; [Romanova et al.|2008| hereafter RKL08)' 
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The most striking result is the fact that in the unstable regime, 
an ordered structure of two tongues forms, and rotates with 
the angular velocity of the inner disk, showing clear high- 
frequency QPO peaks associated with this rotatiorQ Unsta- 
ble accretion with a preferred number of tongues as a pos- 
sible origin for QPOs has been discussed earlier by |Li &| 
|Narayan| ( |2004| ) . In the present research we have shown that 
an even more ordered situation is possible with small mag- 
netospheres, where two ordered tongues rotate with the fre- 
quency of the inner disk. Variation of the accretion rate leads 
to drifting of this QPO frequency. Correlation between the fre- 
quency and accretion rate has been observed in a number of 
accreting millisecond pulsars (e.g., v an der Klis|2000p . Below 
we describe this and other results in detail. In §2 we discuss 
a possible classification scheme for magnetic boundary lay- 
ers. In §3 we describe our model and reference parameters. 
In §4 we consider MBLs during stable accretion. In §5 we de- 
scribe MBLs and formation of QPOs in the unstable regime. 
We present some discussion and conclusions in §6. 



2 CLASSIFICATION OF BOUNDARY LAYERS 

In accreting magnetized stars the disk is stopped by the stellar 
magnetosphere at the magnetospheric, or truncation, radius 
(e.g., |Elsner & Lamb|1977P : 

where M, and fi are the mass and magnetic moment of the 
star, M is the accretion rate through the disk, and the coeffi- 
cient kA is of the order of unity (e.g.,|Long et aH 2005||Besso-| 
|laz et al. 2008). The magnetospheric radius may be small if 
the magnetic moment of the star fi is small, or if the accre- 
tion rate is high. At a given magnetic moment of the star, 
variation of the accretion rate will lead to expansion or com- 
pression of the magnetosphere, and so at very high accretion 
rates the magnetosphere may be completely buried by accret- 
ing matter (e.g., [Gumming et al.|2001[|Lovelace et al.|2005p , 
while at very low accretion rates it strongly expands, and dif- 
ferent types of accretion are expected at different M (e.g., 
[Romanova et al. 20 08 ). An important parameter of the prob- 
lem is the size of the magnetosphere in units of the stellar 
radius, rm/R*- Depending on this parameter we can define 
three t5^es of boundary layers: 

(1). Hydrodynamic BL. If the ratio rm/R* « 1, then the 
magnetic field does not influence the dynamics of matter 
flow in the vicinity of the star, and a purely hydrodjmamic 
boundary layer is expected where the magnetic field can be 
neglected. In this case the disk matter comes to the surface 
of the star and interacts with the star in the equatorial zone 
(e.g., iLynden-B ell & Pringlel|1974| |Pringle|[T9tT) [Popham &| 
[Nara yan 1995 ). Such interaction leads to release of a signifi- 
cant amount of energy at the surface of the star, about a half 
of the gravitational energy, Lbl ~ 0.5GA/,A//7?». Interac- 
tion through such a hydrodynamic boundary layer is expected 
in many non-magnetic white dwarfs and neutron stars and in 
cases where accretion rate is so high that the magnetic field 
is buried by the accreting matter. It has been noted that the 
boundary layer matter will not interact with the star as a thin 



layer, but will spread to higher latitudes along the surface of 
the star due to matter pressu re (jPerland et al.||1982 . This 
effect has been calculated by |Inogamov & Sunyaev| I 1999| ) 
and Piro & Bildsten (2004) and has recently been observed 
in axisjmmetric hydrodynamic simulations by jFisker & Bal-| 
|sara| ( |2005[ ). Much less work has been done for analysis of 
the magnetic boundary layers (MBL) . 

(2) . Magnetic Boundary Layer (MBL) With Magnetospheric 
Gap. If the magnetosphere is only slightly larger than the ra- 
dius of the star, rm/R* > 1, then it stops the disk at a small 
distance from the star. The simulations described in this pa- 
per show that in this case, accretion may be either in the 
stable regime producing funnel streams, or in the unstable 
regime where most of the matter accretes through instabilities 
(KR08a; RKL08). We find that an unusual type of MBL forms 
in the unstable regime with two ordered streams rotating in 
the equatorial plane at approximately the angular velocity of 
the disk. These streams may be responsible for high-frequency 
QPOs. 

(3) . MBL Without a Magnetospheric Gap. If the magneto- 
spheric radius is comparable to the stellar radius, Vm/ R* ~ 1, 
then the magnetosphere cannot stop the disk, and the disk 
interacts with the stellar surface. In this case the magnetic 
field may be somewhat enhanced in the accretion disk due to 
wrapping of the magnetic field lines by the disk matter. Pre- 
liminary simulations have shown that interesting phenomena 
may happen at the disk-magnetosphere boundary, such as in- 
teraction through a Kelvin-Helmholtz-type instability. 

Global 3D simulations of accretion to stars with large 
magnetospheres (a few stellar radii in size) have been per- 
formed earlier (e.g., RUKL03,04; KR05). It was recently found 
that stars may be either in the stable or unstable regime of 
accretion (KR08a,b; RKL08). In this paper we show results 
of global 3D simulations of accretion to star with very small 
magnetospheres, in both stable and unstable regimes, and 
some results for direct interaction between the disk and the 
star. 



3 MODEL 

We solve the 3D MHD equations presented, e.g., in RUKL03, 
using the "cubed sphere" second order Godunov-tj^pe code 
described in |Koldoba et al.| ( |2002p . The equations are written 
in a coordinate system rotating with a star. The energy equa- 
tion is solved for the entropy, and the equation of state is that 
for an ideal gas. Viscosity is incorporated into the code with 
a viscosity coefficient proportional to a l [Shakura & Sunyaev] 
|1973), where a < 1. The action of the viscosity is restricted to 
regions of high density (the disk) . Viscosity helps bring matter 
towards the star from the disk. 

3.1 Dimensionalization and Reference Values 

The MHD equations are solved in dimensionless form so that 
the results can be readily applied to different types of stars. 
We have the freedom to choose three parameters from which 
the reference scales are derived, and we choose those to 
be the star's mass M, radius 7?, and magnetic moment /x* 
(with corresponding magnetic field B, = p,t/Rl). We take 
the reference mass Mo to be the mass M of the star. The 



^ See animations at 
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reference radius is about three times the radius of the star, 
R() = i?./0.35 « 3 X 7?.. In our past research (RUKL03; 
RUKL04) this radius approximately corresponded to the trun- 
cation (or magnetospheric) radius. In the present research 
the truncation radius is much smaller, but we keep this unit 
for consistency with our prior work. The reference velocity 
is vo = (GM/Ro)^^^. The reference time-scale to = Ro/yo, 
and the reference angular velocity coo = 1/to- We measure 
time in units of Po — 2-Kto (which is the Keplerian rotation 
period at r = Rq). In the plots we use the dimensionless time 
T = t/Po. 

The reference magnetic moment = fJ-t/Ji, where Jl is 
a dimensionless parameter which determines the dimension- 
less size of the magnetosphere. The reference magnetic field 
is Bo = fio/Ro- The reference density is po = Bq/vq. The ref- 
erence accretion rate is Mq = po^^o^o- Taking into account re- 
lationships for po. Bo, Vo we obtain Mo = A/Qj^voRo)- The 

dimensional accretion rate is M = MMo, where M is the di- 
mensionless accretion rate onto the surface of the star which 
is obtained from the simulations. Substituting parameters of 
the star to the reference values, we obtain the dimensional 
accretion rate: 



yP J (GAf)i/2(i?,/0.35)7/2- ^ ^ 

One can see that at fixed parameters of the star (Af, Rt, /i,), 
the accretion rate depends on the ratio M/Jl^. Further analy- 
sis shows that parameter Jl^ varies more strongly than M and 
determines the variation of the accretion rate. In the paper we 
vary parameter p which leads to variation of the dimensional 
accretion rate. Table [l] shows examples of reference variables 
for different stars. We solve the MHD equations using the nor- 
malized variables p = p/po, v = v/vo, B — B/Bo, etc. Most 
of the plots show the normalized variables (with the tildes 
implicit). To obtain dimensional values one needs to multiply 
values from the plots by the corresponding reference values 
from Table [l] 



3.2 Initial and Boundary Conditions 

Here we briefly summarize the initial and boundary con- 
ditions, described in detail in our earlier papers (e.g. in 
RUKL03). The simulation region consists of the accretion disk, 
the corona and the star. The star has a dipole magnetic field 
which is frozen into its surface. The disk is cold and dense. 
The corona is 100 times hotter and 100 times less dense (at 
the fiducial point at the inner edge of the disk). Initially we 
rotate the disk with Keplerian velocity and the star with the 
angular velocity of the inner disk, so as to avoid discontinuity 
of the magnetic field lines at the disk-star boundary. We also 
rotate the corona above the disk with the Keplerian veloc- 
ity of the disk so as to avoid discontinuity at the disk-corona 
boundary. This is necessary since the discontinuity may lead 
to artificially strong forces on the disk and strong magnetic 
braking. In addition, we derive such a distribution of the den- 
sity and pressure in the disk and corona, that the sum of the 
gravitational, pressure and centrifugal forces is equal to zero 
ever5where in the simulation region ( [Romanova et al.|2002p . 
This initial condition ensures a very smooth gradual start-up, 
where the disk matter accretes slowly inwards due to viscous 



stresses. On the other hand, this condition dictates a density 
distribution that does not correspond to that in an equilibrium 
viscous disk (e.g. |Shakura & Sunyaev|1973| l. That is why our 
disk starts reconstructing itself on the viscous time-scale. For- 
tunately, matter from the inner parts of the disk flows inward 
for a long time, and the matter flux increases with a, though 
dependence is not linear (as one would expect from the theo- 
retical prediction for viscous disks) . In addition, interaction of 
the disk with the external magnetosphere leads to some mag- 
netic braking, leading to an enhancement of the accretion rate 
which may be smaller or larger than the effect of the viscos- 
ity and also to oscillations of the matter flux. That is why we 

cannot use the dimensionless accretion rate M as a parameter 
of the problem. Instead, we choose p as the main parameter 
which is responsible for the size of the magnetosphere. Equa- 
tion (1) shows that the dimensional accretion rate onto the 

star is regulated by the combination M /p^. In this paper we 
consider small magnetospheres with p ~ 0.08 — 0.5, which 
is about 10 times smaller than in our past research (RUKL03; 
RUKL04; KR05). 

The size of the simulation region is Rmax = 15-Ro ~ 
45-R, . We initially place the inner radius of the disk at a dis- 
tance rd from the star which is either 1.2_Ro ~ 3.67?, (in sta- 
ble cases) or 1.8-Ro ^ 4.87?, (in the unstable cases). Simu- 
lations show that the result does not depend on ra, because 
in both cases the disk matter moves inward and is stopped 
by the magnetosphere at small distances from the star. Then 
we have the freedom to spin the star up or down gradually. 
In this paper we choose a slowly rotating star with angular 
velocity lj, = (GM jrcorY^'^ corresponding to the corotation 
radius Vcor = 2i?o. In dimensionless units, uJl = 0.354, and in 
application to neutron stars, P, = 6.2ms. 

The boundary conditions are similar to those used in 
our other papers (e.g., RUKL03). On the star (the inner 
boundary) the magnetic field is frozen into the surface of the 
star (the normal component of the field, _B„, is fixed), though 
all three magnetic field components may vary. For all other 
variables A, free conditions are prescribed: dA/dn — 0. The 
total velocity vector is adjusted to be parallel to the total mag- 
netic field to enhance the frozen-in condition. At the external 
boundary, r — Rout, again free boundary conditions are pre- 
scribed for all variables. In some cases we fix the density on 
the external boundary so that matter of the disk stays in the 
region during the viscous reconstruction mentioned above. 
We do not supply matter to the disk from the external bound- 
ary, because the amount of matter in the inner regions of the 
disk is sufficient to support accretion over the duration of our 
simulations. 

The grid. Simulations of accretion to stars with small 
magnetospheres require a finer grid compared to the larger 
magnetosphere cases (e.g. RUKL03). We use a grid resolution 
of iVr X = 100 X 41^ (in each of the 6 blocks of the cubed 
sphere) for most of our cases, and A'^,. x A^^ = 130 x 61^ 
for the smallest magnetospheres. The simulations were done 
using 60-180 processors each on the NASA HPC clusters. 

4 ACCRETION IN THE STABLE REGIME 

Below we describe the results of simulations of stable accre- 
tion to stars with small magnetospheres p = 0.08 — 0.2. 
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CTTSs White dwarfs Neutron stars 



M(Mq) 0.8 1 1.4 

R, 2Rq 5000 km 10 km 

(G) 103 106 3 X 10^ 

Ro (cm) 4 X 10" 1.4 x 10^ 2.9 x 10^ 

uo (cm s-i) 1.6 X lO'^ 3 X 10** 8.1 x 10^ 

a;o (s-i) 4 X 10-5 0.21 2.8 X 10^ 

vo 0.55 day-i 3.2 x lO^^ Hz 4.5 x 10^ Hz 

Po 1.8 days 29 s 2.2 ms 



Table 1. Sample reference values of the dynamical quantities for different types of stars. We choose the mass, radius and magnetic field of the 
star, and define the other variables in terms of these three quantities. Note that the frequency j/q is the Keplerian frequency at radius Rq. 




X X 



Figure 1. Stable accretion to a star with a small magnetosphere, 
where the magnetic moment of the star is ^ = 0.2 and is misaligned 
relative to the rotational axis f2 at 6 = 15°. Leftpanel: Density distri- 
bution and sample magnetic field lines in the xz ifi — Q) plane. Right 
panel: The same, but in the xy— (equatorial) plane. The dimension- 
less density varies between 0.8 (red) and 0.005 (blue). 



Simulations show that the boundary between the sta- 
ble and unstable regimes depends on a number of factors, 
such as the accretion rate (determined by the viscosity pa- 
rameter a), the rotation period of the star P,, and the mis- 
alignment angle 6 of the dipole (KR08a,b). If a star with a 
small misalignment angle, <d < 5°, rotates slowly (the disk 
rotates much faster than the star), then the accretion has a 
tendency to be unstable irrespective of a (KR08a,b). We have 
chosen slowly rotating stars with dimensionless angular ve- 
locity widetildeujt = 0.354 (corresponding to a corotation 
radius of rcor ~ 2) . We take slowly rotating stars in order to 
model the situation in which the star initially has a large mag- 
netosphere and is in the rotational equilibrium state, but later 
the accretion rate increases and the magnetosphere becomes 
small. To stabilize accretion we chose & — 15°. We also used 
a = 0.04. In the description of the results we measure time 
in units of the Keplerian rotation period at r = 1. We con- 
sider two main cases, one with Jl — 0.2 which has a small 
magnetosphere, and the other with ^ = 0.08 which has a tiny 
magnetosphere. 

Fig. [l] shows an example of accretion to a star with 
Jl — 0.2. A small magnetosphere forms, in which the mag- 
netic field of the dipole dominates. The position of the inner 
disk edge is determined by the balance between the magnetic 
and kinetic matter pressure. Matter is lifted above the mag- 
netosphere forming two funnel streams, and accretes to the 
surface of the star forming two hot spots. The hot spots show 
the distribution of the total energy flux on the surface of the 
star (see RUKL04 and KR05 for details). 



The hot spots have a preferred position, although in case 
of small misalignment angles the funnel streams are dragged 
by the rapidly rotating disk and may rotate faster than the 
star. Such spot rotation has been observed in earlier simula- 
tions (RUKL03; [Romanova et al.|[2006| . In the opposite sit- 
uation when the star spins fast, both funnels and spots may 
rotate slower than the star (e.g. Romanova et al. 2002). Faster 
or slower rotation of spots may lead to a QPO feature with fre- 
quency higher or lower than the stellar frequency (Romanova 
et al. 2006). It has been suggested that at small misalignment 
angles, the spots may wander around the magnetic poles, pos- 
sibly causing intermittency in some millisecond pulsars Lamb 
et al. (2006 lb) (see also (Altamir ano et aLpOOSj [Casella 
,et al. 2008 )). On the other hand, if the misalignment angle 
of the dipole is not very small, then such faster or slower spot 
motions become less significant. The spots acquire a preferred 
place on the surface of the star and corotate with the star (see 
also RUKL03; RUKL04; KR05). In simulations with 9 = 15° 
we observe a mixture of spots: on the one hand, the disk drags 
the funnel around the weak magnetosphere in spite of the 
relatively high misalignment angle, and frame to frame anal- 
ysis shows that the spots often move faster than the star. On 
the other hand, the hot spots spend a somewhat longer time 
in the /i — f2 plane compared with other positions, due to 
which we observe a definite peak associated with star's spin 
in the Fourier spectrum. Longer simulation runs are probably 
required to see the QPO peak. 

Fig. [2] shows an example of accretion to a star with a 
tiny magnetosphere, ^ = 0.08. It is amazing that in this case 
too, the magnetosphere channels the accreting matter, form- 
ing tiny funnel streams hitting the surface of the star. Density 
waves form in the equatorial plane of the disk. Sometimes 
density fluctuations in the accreting matter push the disk to 
the stellar surface. However, later the magnetosphere is re- 
stored again. The hot spots have a preferred position but they 
often rotate faster than the star as a result of the difference in 
angular velocities between the disk and the star. In this case 
again, a QPO at the stellar frequency is expected, and also a 
QPO corresponding to the frequency at the inner edge of the 
disk and beats between these two frequencies. Longer simu- 
lation runs are required to extract these frequencies. 

The main result of this section is that even small and 
tiny magnetospheres disrupt the disk and channel matter to 
the star, forming hot spots. It is also important to note that 
the high-frequency QPO associated with rotation of the inner 
edge of the disk is expected if the disk rotates faster (slower) 
than the star. Signs of faster rotation of the spot are clearly 
observed. Drifting of the high-frequency QPO is expected if 
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Figure 2. Stable accretion to a star witii a very small (tiny) magnetosphere (p = 0.08) with = 15° (o = 0.04). Top two panels: Density 
distribution and sample magnetic field lines in the xz {fi — H) plane and the xy— (equatorial) plane in a coordinate system rotating with the star. 
The dimensionless density varies between 1.4 (red) and 0.006 (blue). Bottom panel: Distribution of the dimensionless energy flux of matter onto 
the surface of the star (pole-on), varying between 0.4 (red) and 0.0009 (blue). The time T is measured in units of AT = I/I6P0, where Po is 
the Keplerian period at r = 1. 



the inner edge of the disk varies as a result of variation of the 
accretion rate. Future (longer) simulations will help obtain 
different frequencies. 



5 ACCRETION IN THE UNSTABLE REGIME: A NEW TYPE 
OF BOUNDARY LAYER AND QPOS 

5.1 Formation of two symmetric streams and spots 

Now we take a star with a small misalignment angle, 6 = 5°, 
and investigate accretion in the unstable regime. First we 
choose a small magnetosphere with magnetic moment Jl = 
0.2. We take the viscosity coefficient a — 0.1 in all runs be- 
low. A star rotates slowly with angular velocity uu — 0.354 
(corresponding to rear = 2). We observe that the instability 
starts, and matter penetrates through the equatorial region of 
the magnetosphere through a number of tongues. Later, how- 
ever, two sjTTimetric ordered tongues form and rotate syn- 
chronously with angular velocity approximately equal to the 
angular velocity in the inner region of the accretion disk, that 
is, much faster than the star. These tongues reach the surface 
of the star and produce two hot spots at the star's equator 



which move faster than the star. The tongues deposit a signif- 
icant amount of energy to the surface of the star Part of this is 
the gravitational energy associated with acceleration of mat- 
ter towards the star. We take this energy into account while 
plotting hot spots. In addition, there is energy released due 
to friction between the surface of the slowly rotating star and 
the foot-points of the rapidly rotating tongue. 

Fig. [3] shows snapshots of rotation of the tongues, slices 
in the equatorial plane and hot spots (energy flux) on the 
star's surface. The slices show that the magnetosphere has a 
strongly modified shape, particularly in those places where 
the tongues reach the surface of the star. They push the mag- 
netosphere to the surface of the star and when the tongue 
moves, the magnetosphere re-emerges, while the tongue 
pushes another part of the magnetosphere to the surface of 
the star Two strong hot spots form close to equatorial region. 
Sometimes they are not symmetric because the tongues push 
the magnetosphere more to one side than another Some mat- 
ter accretes to the poles through weak funnel streams. It is 
interesting that the polar spots also rotate with the angular 
velocity of the tongues. Thus, in the case of a small magneto- 
sphere we observe a new phenomenon - a modified boundaiy 
layer, where the funnel streams and hot spots move with the 
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Figure 3. Snapshots of MBL accretion througii the synchronized instabiUty tongues in the case of a relatively large magnetosphere (p, = 0.2 in 
dimensionless units). The time interval shown is a small portion of the simulation; the tongue pattern is steady for the entire duration of the 
simulation. Top panel: A surface of constant density (0.4 in dimensionless units) and sample magnetic field lines. Middle panel: Density distribution 
and sample magnetic field lines in the ay-plane. The density varies from 6.7 (red) to 0.01 (blue). Bottom panel: Energy flux onto the star's surface, 
ranging from 3.4 (red) to 0.008 (blue). The figures are shown in a coordinate system rotating with the star The time T is measured in periods of 
Keplerian rotation at r = 1. 



angular velocity of the disk — much faster than the star. This 
is a new type of boundary layer where matter of the disk in- 
teracts with the star through unusual symmetric tongues. Ro- 
tation of the hot spots along the surface of the star is expected 
to give strong QPO peaks at twice the frequency of the inner 
disk. 

To check this phenomenon, we performed a set of simula- 
tion runs at a variety of stellar magnetic moments: Jl — 0.16, 
]! = 0.15, Jl — 0.12, Jl = 0.1, keeping all the other parame- 
ters fixed. We observed that similar symmetric streams form 
and rotate around the star. However, at smaller Ji the disk 
stops closer to the star. In the case of fl — 0.1, only a tiny 
magnetosphere forms with tiny streams (see Fig. |4| . The bot- 
tom panels of the Fig. |4] show the density distribution in the 
hot spots instead of the emitted flux, because it is expected 
that energy would be released mainly due to friction between 
the tongues and the surface of the star, which we do not cal- 
culate in this paper. In spite of the fact that the tongues are 
very small, there is still energy associated with gravitational 
acceleration of matter in the tongues and the corresponding 
heating of the stellar surface. The hot spots associated with 
gravitational acceleration are located in the regions of high- 



est density and occupy a much smaller area than the spots 
shown in Fig.[4| At the end of the ^ = 0.1 simulation run the 
accretion rate increased, the disk reached the surface of the 
star and started interacting with the star's surface through an 
equatorial belt (see §4.3). 

We also performed exploratory simulations of accretion 
to stars with larger magnetospheres, Jl = 0.5 and observed 
that in many cases accretion through two or one ordered 
tongues dominates. We varied the a parameter (which reg- 
ulates the accretion rate) and rotation rate of the star and 
obtained tongues and hot spots with different levels of order. 
In some cases two identical, diametrically opposite streams 
dominate accretion during the whole simulation run. In other 
cases only one stream dominates while the other is weaker, 
and therefore only one hot spot will determine the frequency 
observed in the light curve. In yet another type of case, mul- 
tiple tongues are observed (like in cases of large magneto- 
spheres, e.g. KR08a; KR08b), but one or two tongues are 
stronger than others, and can give rise to a QPO peak. Be- 
low we include the Jl = 0.5 case as an example of accre- 
tion to slightly larger magnetospheres. This case links the very 
small magnetospheres considered in this paper with the much 
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T=3.38" 




Figure 4. Top panel: A surface of constant density (0.47 in dimensionless units) and sample magnetic field lines in the case of a small magneto- 
sphere, Jl = 0.1. Middle panel: xy-slice showing the density distribution and sample magnetic field lines. Bottom panel: density distribution on 
the surface of the star, ranging from 3.8 (red) to 0.01 (blue). The figures are shown in a coordinate system rotating with the star. The time T is 
measured in periods of Keplerian rotation at r = 1. 



larger ones considered in KR08b, where QPOs of similar ori- 
gin are observed at some sets of parameters. 



5.2 Frequency analysis 

We performed frequency analysis for all the above cases. First 
we perform what we call spot-omega analysis, which we find 
to be the most informative for analysis of our simulations. 
Namely, we plot the equatorial distribution of the emitted en- 
ergy flux at different moments of time (see Fig. [sj left col- 
umn). We obtain diagonal lines which show that the spots 
have a definite order in their rotation along the star's equator 
The slopes of these lines are proportional to the angular veloc- 
ity of the spots. We obtained a number of almost parallel lines 
which reflect multiple rotations of a single spot with approx- 
imately the same angular velocity. We performed this spot- 
omega analysis for all cases. Fig.[5]shows the spot-omega dia- 
grams for the same time interval. One can see that at smaller 
fl, the lines are steeper because the inner disk is closer to the 
star and the rotation of the streams/spots is faster Near the 
lines we show the rotation frequencies of the corresponding 
spots in units of the stellar rotation frequency. Note that in the 
case of the smallest magnetosphere (Jl = 0.1), the streams 
become very small and the spots associated with the energy 
flux due to gravitational acceleration become weak. In this 
case we plot the distribution of the density along the equator 
Closer to the end of this simulation run, the disk comes to the 
surface of the star, and the streams (and ordered lines) disap- 
pear (see bottom left panel of Fig.[5]l. In the case of the largest 



magnetosphere, fl — 0.5, the spot-omega diagram shows a lit- 
tle less order because the magnetosphere is stronger and the 
rotation of the tongues is less synchronized. 

Next we performed wavelet analysis of the light-curves 
from hot spots. In the cases with Jl = 0.5, 0.2, 0.15 and 0.12, 
the light-curve is calculated with the suggestion that all the 
gravitational and thermal energy of the matter falling onto 
the star is converted into thermal energy, which is radiated 
isotropically. The light from the surface of the star is inte- 
grated in the direction of the observer, whose line of sight 
makes an angle of i = 90° with the rotational axis of the star 
(see RUKL04 and KR05 for details). In the case with Jl ^ 0.1, 
it is expected that energy is released mainly due to friction be- 
tween the rapidly moving streams and the star Calculation of 
the radiation from such friction is beyond of the scope of this 
paper As a first approximation we suggest that the emitted 
flux is proportional to the matter density at the star's surface. 
Fig. [s] shows wavelet spectra for all these cases. These plots 
exclude early times, when the streams have yet to reach the 
star's surface. Since the wavelet transform uses a window of 
width At centered at time t to calculate the time-dependent 
frequency spectrum of the lightcurve from time ti to time t2, 
it is necessary to have ti + At/2 < t < t2 — At/2. This results 
in a portion of the wavelet plot being cut off. The width At 
is inversely related to the frequency, so lower frequencies are 
more strongly affected by this restriction. One can see that in 
each case several frequencies are observed. Comparisons of 
the wavelet frequencies with the spot-omega frequencies (left 
panels) show that one of the frequencies is approximately 
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Figure 5. Frequency analysis for runs with different magnetospheric sizes. Left column: Spot-omega analysis showing the emitted energy flux 
distribution in the equatorial plane at different times. The flux varies from 3.4 (red) to 0.008 (blue). The lines reflect the motion of individual 
spots on the stellar surface. The slope of the lines is proportional to the angular velocity of the hot spots. Right two columns: Wavelet and Fourier 
spectra of the light curves from the hot spots obtained at an observer inclination angle of j = 90° (i.e., when the observer is in the equatorial 
plane). The arrows show the frequency corresponding to the presence of two rotating spots. In the bottom row the lines show the distribution 
of the density, and not the flux, in the spots, and the density varies between 3.8 (red) and 0.01 (blue). In this case the lightcurve for which the 
wavelet and Fourier spectra are shown is calculated assuming that the emission of the star is proportional to the density. 



twice as high as the frequency of the single spots shown in 
the spot-omega diagram, and hence corresponds to the mo- 
tion of two spots along the surface of the star 

We also performed Fourier analysis of the light-curves. 
For this analysis we chose only the time intervals during 
which steady rotation of the spots was observed. In the case 
with Ji = 0.1, we also excluded the late times after the spots 



disappear and the boundary layer forms. The Fourier spectra 
show similar peaks as the wavelet. Between all these plots, 
the peak corresponding to rotation of two spots is most im- 
portant. Thus we get an approximate agreement between the 
high-frequency QPOs obtained in all three methods of the fre- 
quency analysis. 

There are also lower frequency peaks observed in both 
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the wavelet and Fourier plots. These frequencies may either 
reflect the shape of the spots, or be beat frequencies (see also 
[Smith et al.|1995|[Psaltis et al.|1998| l. We do not see the star's 
rotation frequency itself. However, it is possible that our simu- 
lations are not long enough for the wavelet analysis to capture 
possible QPOs associated with the slowly rotating star On the 
other hand, the spot-omega diagrams do not show signs of 
stellar rotation, so it is possible that the corresponding QPO 
is very weak. 

Note that spiral density waves often form in the disk (see 
Fig. [4] middle panels). The density waves tend to be station- 
ary in the coordinate system rotating with a star, that is, they 
are generated by the inclination of the dipole. The rotating 
tongues disrupt the inner part of this pattern, but the exter- 
nal pattern approximately corotates with the star (see Fig.[4|. 
The figure also shows that sample magnetic field lines start- 
ing out equidistantly from the star's surface also accumulate 
in the density waves. Note that these density waves are simi- 
lar in shape to those suggested by Miller et al. ( 1998 ), except 
that here the streams are guided by magnetic fields instead of 
being produced by radiation drag. 




Figure 7. Enlarged view of two snapshots from Fig. |6] showing a 
surface of constant density (0.43 in dimensionless units) and sample 
magnetic field lines. The top and bottom panels show different pro- 
jections. One can see that matter is lifted along the surface of the star 
to larger latitudes. In the case when the disk comes to the surface of 
the star, a new instability appears at the disk-star boundary, which is 
probably of the Kelvin-Helmholtz type. 



5.3 Interaction through the boundary layer 

Closer to the end of the simulation run with the smallest mag- 
netosphere (Jl — 0.1), the accretion rate increased and the 
disk matter started interacting with the surface of the star 
along a belt-like path (see Fig. [6]l. Later, the belt became 
wider, and an unusual pattern formed, connected with some 
instability. This is probably the Kelvin-Helmholtz instability 
which develops between the slowly rotating stellar surface 
and the much more rapidly rotating disk, which is expected 
to rotate at a frequency lydisk ~ Si^t at the surface of the star. 

We should note that in all cases, when the modified 
tongue or a boundary layer reaches the surface of the star, 
it spreads to higher latitudes as it was predicted by Inogamov 
& Sunyae v] ([1999P (see also Piro & Bildsten 2004; Fi sker &| 
Balsara,2005 I. Fig. [7] shows such spreading and also a closer 
view of the unstable boundary layer on the star's surface. Dur- 
ing such close contact between the disk and the star, energy is 
expected to be released mainly from friction between the disk 
matter and the stellar surface, like in the usual hydrodjmamic 
boundary layer. Compared to all the cases in the previous sec- 
tions, no energy is associated with the matter falling onto the 
star's surface. Such boundary layers require special investiga- 
tion. Note that even in the case with no magnetosphere, a 
weak magnetic field threads the disk, and magnetic field lines 
are wrapped in the inner parts of the disk. 

At later times, when the disk comes closer to the star, 
the magnetic field lines trapped in the equatorial region are 
pushed closer to the star or buried. Note that at the same time 
the field lines above and below the disk are not buried (see 
also Fig. [7]l. Possibly the process of field burial by the thin 
disk through the boundary later should be reconsidered tak- 
ing into account the fact that a significant amount of magnetic 
flux may inflate into the corona and stay there. 

5.4 Example: Application to Accreting Millisecond 
Pulsars 

Here we show a sample application of our simulations to ac- 
creting millisecond pulsars. We take a neutron star with mass 



M* — 1.4M0, radius 7?, = 10 km and surface magnetic 
field B, = 3 X 10* G. The corresponding reference values 
are in Table [l] We convert the dimensionless results obtained 
in §4.1 and §4.2 into dimensional values. The dimensionless 
stellar angular velocity is = 0.354. The corresponding di- 
mensional angular velocity is uj, = uuujo — 1002 s"\ di- 
mensional period is P» = 27r/a;* = 6.3 ms, and dimensional 
frequency is — 159 Hz. The time T used in the figures 
has the reference value Po — 2.2 ms which is the Keplerian 
rotation period at r = 1 (3i?* = 30 km). The spot-omega 
diagrams in Fig. [s] (left panels) show the frequencies associ- 
ated with the rotation of one spot for different values of Jl. In 
the case of the larger magnetosphere (fl — 0.2) the frequency 
is uispot ~ 3Aiyt — 541 Hz. Since there are two antipodal 
hot spots, we expect the observer to see twice that frequency, 
i^2spots ~ 1082 Hz. In the case with Jl — 0.15, the disk is 
closer to the star and the frequencies are uispot ~ 4.5i^* ^ 715 
Hz and iy2spots ~ 1430 Hz. For Jl = 0.12 the frequencies are 
J^ispot ~ 5.3i/* — 843 Hz and U2spots ~ 1686 Hz. And for 
Jl — 0.1 the frequencies are vispot ~ 6.2i^, — 986 Hz and 
V2spots ~ 1972 Hz. The frequencies are quite high because 
the disk is close to the surface of the star In simulations with 
a larger magnetosphere, Ji = 0.5, we obtain lower frequen- 
cies, uispot ~ 2.2i^, = 350 Hz and i'2spots ~ 700 Hz. 

In application to neutron stars the formulae (1) can be 
re-written as: 



M = 2.2 X 10" 



2.8 X 1033g 



B, 



3 X lO^G 



106cm 



Mq 



(2) 



One can see that the accretion rate depends on the dimension- 



less parameter M //i^. Table[2]shows that the ratio M /Jl'' de- 
creases systematically with Jl, and so do the accretion rate and 
the frequencies vispot and V2spots- So, the considered above 
cases with different Jl correspond to different accretion rates 
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Figure 6. Same plot as in Fig. [4] but closer to the end of the run when the disk approaches the star. The surface shown in the top panels has a 
density of 0.43. In the bottom panels the density varies between 2.6 (red) and 0.007 (blue). The figures are shown in a coordinate system rotating 
with the star. The time T is measured in periods of Keplerian rotation at r = 1 . 





/I = 0.5 


Jl = 0.2 


Jl = 0.16 


Jl = 0.15 


Jl = 0.12 


Jl = 0.1 




0.88 


3.0 


4.1 


4.4 


5.5 


6 


M (Mgyr-i) 


1.9 X 10"^ 


6.6 X lO"'' 


9.0 X 10-8 


9.7 X 10"'' 


1.2 X 10"* 


1.3 X 10"** 


l^lspot (Hz) 


350 


541 


670 


715 


843 


986 


I'2spot (Hz) 


700 


1082 


1340 


1430 


1686 


1972 



Table 2. Values of the quantities described in §5.4 as a function of the magnetospheric size parameter /j.. 



where, depending on accretion rate, the frequency drifts be- 
tween V2spots ~ 700 Hz when the accretion rate is lower and 
the magnetosphere is larger, and v-zspota ~ 1972 Hz when the 
accretion rate is higher and the magnetosphere is very small. 
The ratio between these two frequencies is 2.8 and can be 
smaller or possibly larger depending on the variation of the 
accretion rate, and is close to the observed drifts of the main 
high-frequency QPO in, e.g., millisecond pulsars {van der Klis| 
|2000| ) and dwarf novae ( [Warner & Woudt|2006| ). 



Only if the two hotspots are antipodal and identical, do 
we expect uispot to be absent from the power spectrum. This 
is true for the small magnetosphere (Jl < 0.2) cases. In the 
large magnetosphere (Jl = 0.5) case, we see that the spots 
are not identical; one spot is often much larger than the other, 
and hence the frequency uispot may dominate. If the magnetic 
field of the star is not an ideal dipole field (that is, if it is a 
slightly misplaced dipole, or a more complex field (e.g., |Long| 
|et aL|20()7||2008 ) then the symmetry breaks and one spot will 
be always larger than the others, and the lower frequency 



i^ispot will dominate. Then we expected the frequency frift 
between 350 Hz and 990 Hz. 

Simulations of stars with larger magnetosphere sizes, 
Jl > 0.5, usually show much more chaotic behavior in the 
unstable regime (KR08a,b; RKL08). In spite of stochastic ac- 
cretion in the unstable regime, the spots on the surface of 
the star show a component which rotates with the angular 
velocity of the inner disk (KR08b) . This component analyzed 
through rotation of spots (like in the left column of Fig. [sj 
shows clear inclination of lines associated with motion of dif- 
ferent temporary spots with angular velocity approximately 
equal to the angular velocity of the inner disk (KR08a) . This 
frequency component is weaker than the stochastic compo- 
nents associated with unstable accretion. However in longer 
simulation runs this QPO peak associated with the inner disk 
frequency may by amplified with time and may become signif- 
icant, because the high-amplitude stochastic components are 
relatively incoherent and their frequencies constantly drift. 

Although we consider accretion to a very slowly rotating 
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star, similar MBLs and QPOs are expected for more rapidly 
rotating stars. 



5.5 Disk oscillations 

Note that some peaks observed in the wavelets and Fourier 
spectra may be associated with the disk-star interaction. An 
accretion disk can have different modes of oscillation. These 
include bending oscillations of the inner disk driven by the 
star's rotating misaligned dipole field (e.g., Lai & Zhang 2008 ) 
and radial oscillations of the inne r disk (e.g., ,Alpar & PsaltiS| 
2008: Lov elace & Romanova|2007t|Lovelace et al.|2009||Erkut| 
et al.||2008 l. The bending oscillations lead to the formation 
of an m = 1 mode spiral wave rotating with the frequency 
of the star. In our simulations we do see traces of bending 
waves. The radial oscillations can arise from a linear Rossby 
wave instability (RWI) where the angular frequency of the 
oscillations (with mode number m = 1) is less than the peak 
in the angular velocity in the inner disk, Slmax (Lovelace &[ 
|Romanova|2007[ [Lovelace et al.|2009P . This mode is radially 
trapped in a narrow region inside the radius at which f2max 
occurs. The beat between this mode and the stellar rotation 
frequency may lead to the coupled twin peak QPOs observed 
in some millisecond pulsars (e.g., van de r Klis 2006 ). 

Here we derive from our simulations the radius at which 
the disk angular velocity matches that of the spot. Fig. [8] 
shows this analysis for cases with different Jl. First, we plot 
the radial dependence of the angular velocity in the disk at 
sample times (see the top panels of the Fig. [8]l. Then we 
take the angular velocity of the spot from the spot-omega 
diagram (Fig. [s] left panels) and find the point where an- 
gular velocity of the spot equals the angular velocity in the 
disk (see red dots). Next, we plot the radial dependence of 
angular momentum (Fig. [8] middle panels) and density (bot- 
tom panels) at same moment of time. The circular line in the 
plots is the line on which the angular velocity equals that 
of the spot. One can see that in all cases the rotational ve- 
locity of the spot corresponds to the inner edge of the disk. 
It corresponds to a distance only slightly larger than the 
Alfven surface (red line), where the kinetic plasma parameter 
l3i = {p + pv'^)/{B'^/8tv) = 1. Note that the magnetosphere is 
strongly non-axisjmimetric, and therefore the distribution 
in the disk varies with time. 

Different oscillation modes in the disk may influence the 
position or brightness variation of the spots on the surface of 
the star. We do observe some frequencies in the power spec- 
tra of the light-curves, but at present we are not sure that the 
observed frequencies reflect disk oscillations. Future (longer) 
simulation runs may help reveal such a possibility. On the 
other hand, disk oscillations can be investigated directly from 
the simulations. We plan to do this in the future. 

Fig. [5] shows that there are several peaks observed in 
the wavelet and Fourier spectra. We know the origin of only 
one of them, associated with rotation of the unstable tongues. 
Other peaks may be associated with disk oscillations or with 
beat frequencies between the disk and the star. However, at 
present our runs are not long enough to establish the origin 
of these peaks. 



6 CONCLUSIONS 

We considered accretion to slowly rotating stars with small 
and tiny magnetospheres in the stable and unstable regimes. 
We conclude that: 

(1) In the stable regime, a small magnetospheric cavity 
and tiny funnel streams form, producing hot spots on the 
star's surface which tend to be in a preferred position deter- 
mined by the dipole inclination (we had example runs for 
6 = 15°). In this case the frequency of the star is expected to 
dominate. However, the rapidly rotating disk has a tendency 
to "drag" the funnel stream to faster rotation, due to which 
parts of the hot spots often rotate much faster than the star. 
At lower O this effect becomes even more significant, and the 
spots may rotate faster than the star for a long time, leading 
to QPOs (see also RUKL04), or slower than the star leading to 
matter accreting through a trailing funnel, (e.g., Romanova et 
al. 2002) producing a lower-frequency QPO. 

(2) In the unstable regime we observed that matter ac- 
cretes through two ordered streams which rotate with the an- 
gular velocity of the inner disk, that is, much faster than the 
star. They hit the surface of the star forming two antipodal 
hot spots. Rotation of these spots along the surface of the star 
leads to high-frequency QPOs. Such persistent streams/spots 
have been observed at a variety of parameters and are seen 
to be quite long-lived. Coherent rotation of this kind has 
been observed for small magnetospheres. For large magne- 
tospheres we find that the spots are much more chaotic in 
the sense that they form at different parts of the star (RXL08; 
KR08a). In intermediate cases we observed that one or two 
ordered streams may form, which are less ordered than in the 
cases of small magnetospheres, but still may give QPO peaks 
(KR08b). On the other hand, the accretion may be stochastic, 
but one or two streams may be stronger than the others, and 
the hot spots associated with these streams may lead to QPOs. 

(3) Correlation is observed between the size of the mag- 
netosphere and the QPO frequency: the frequency is higher 
at smaller magnetospheric sizes. We expect that secular vari- 
ation of the accretion rate will lead to drifting of the QPO fre- 
quency. Correlation between the frequency and the accretion 
rate has been observed in a number of accreting millisecond 
pulsars (e.g., [van der Klis|2000[ ). 

(4) We were able to model a wide range of QPO frequen- 
cies. In application to millisecond pulsars the frequency asso- 
ciated with rotation of one spot varies between 1^1 spot = 350 
Hz and 990 Hz. In cases of small magnetospheres, the two 
spots are very similar in brightness, and the expected frequen- 
cies are twice as high. Only if the two spots are antipodal and 
identical, do we expect uispot to be absent from the power 
spectrum. This is true for the small magnetosphere QI < 0.2) 
cases. In the large magnetosphere (/i=0.5) case, we expect to 
see lyispot- 

The most striking result of this paper is the discovery of 
a new regime of unstable accretion which shows clear high- 
frequency QPO peaks. 
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Figure 8. Top panels: Radial dependence of tiie disk angular velocity for the moments of time corresponding to the lower panels (thick blue line) 
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thick line) and the circle shows where the angular velocity equals that of the spot (black dash-dot line) .Bottom panels: Same, but for the density 
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